10 REM - FROM: 'FOURIER SMOOTHING WITHOUT THE FAST FOURIER TRANSFORM',
20 REM -   BY ERIC E. AUBANEL & KEITH B. OLDHAM ; 'BTYE', FEB 1985, P. 207
30 REM
40 REM - MODIFIED FOR THE COMMODORE-64 BY J.J. WORMSER, APRIL 1985
100 PRINT"[147]"
110 PRINT"    FOURIER SMOOTHING PROGRAM "
120 PRINT:PRINT" THIS PROGRAM CALCULATES THE OPTIMUM"
130 PRINT" SMOOTH CURVE USING THE FOURIER INVERSE"
140 PRINT" TRANSFORM IN A QUADRATIC DIGITAL"
150 PRINT" FILTER PROCESS. THE PROGRAM ALSO"
160 PRINT" INCLUDES A STRAIGHT-LINE CORRECTION"
170 PRINT" ROUTINE WHICH REDUCES 'END EFFECT'.":PRINT
180 PRINT" YOU SIMPLY ENTER THE DATA AS PROMPTED. "
190 PRINT" THE NUMBER OF TRANSFORM POINTS TO BE"
200 PRINT" KEPT ('E') DETERMINES THE DEGREE OF"
210 PRINT" SMOOTHNESS. THIS NUMBER MUST LIE"
220 PRINT" BETWEEN 1 AND INT((N+1)/2). THE"
222 PRINT" DEGREE OF SMOOTHNESS IS INVERSELY"
224 PRINT" PROPORTIONAL TO THIS VALUE."
226 PRINT
230 PRINT" HIT ANY KEY TO CONTINUE";
240 GETK$:IFK$=""THEN240
250 PRINT"[147]":CLR
260 INPUT" ENTER NUMBER OF DATA POINTS";N
270 REM LEAVING R&I ARRAYS UNDIMENSIONED    LIMITS VALID VALUES OF E TO <=10
280 N2=INT((N+1)/2+1):DIMX(N),X1(N),U(N2),V(N2)
290 FOR I=0TON-1
300 INPUT" ENTER DATAPOINT VALUE";X(I)
310 NEXTI
320 GOSUB450
330 PRINT"[147]":PRINT" SMOOTHED DATA FOR 'E'=";ES;"ARE:":PRINT
340 FORI=0TON-1
350 PRINT" X(";I;")=";X1(I)
360 NEXTI
370 PRINT:PRINT" IF YOU WANT TO TRY A DIFFERENT E,"
380 INPUT" ENTER (1) ELSE ENTER (0)";MO
390 IF MO<1ORMO>1THEN420
400 IF  MO=1THENGOSUB450
410 GOTO330
420 END
430 REM FOURIER ALGORITHM SUBROUTINE        BEGINS AT LINE 60. LINE #'S ARE THE
440 REM SAME AS FOR THE HP VERSION OF THE   SUBROUTINE
450 PI=3.141593
460 PRINT"[147]":PRINT" NUMBER OF TRANSFORM POINTS 'E'";
470 INPUT ES
475 PRINT
480 IFES>INT((N+1)/2)THENPRINT" 'E' TOO LARGE":FORI=1TO1000:NEXT:GOTO460
490 IFES<>INT(ES)ORES<=1THENGOTO460
500 IFES<=QTHEN1260
510 REM
520 IFQ=<>0THEN720
530 REM STRAIGHT-LINE CALCULATION
540 S1=0:S2=0
550 D=INT(N/10)
560 FORJ=0TOD-1
570 S1=S1+X(J)
580 S2=S2+X(N-J-1)
590 NEXTJ
600 X1=S1/D:X2=S2/D
610 M=(X2-X1)/(N-D)
620 B=(X1+X2)/2-M*N/2
630 REM CALCULATE(R0)
640 G=0
650 FORJ=0TON-1
660 X(J)=X(J)-M*J-B
670 G=G+X(J)
680 NEXTJ
690 R(0)=G/N
700 Q=1
710 REM
720 PRINT" WORKING ON CALCULATIONS. PLEASE WAIT..."
730 PRINT
740 J2=INT((N-1)/2)
750 P1=INT(100*LOG(2*J2-1)/LOG(2)+.5)/100
760 FORK=QTOES-1
770 J1=J2
780 S=PI*K*2/N
790 C=COS(S):S=SIN(S)
800 FORJ=1TOJ1
810 L=2*J-1
820 U(J)=X(L)*C+X(L+1)
830 V(J)=X(L)*S
840 NEXTJ
850 S=2*S*C:C=2*C*C-1
860 FORP=1TOP1
870 U(J1+1)=0:V(J1+1)=0
880 J1=INT((J1+1)/2)
890 FORJ=1TOJ1
900 L=2*J-1
910 U=U(L)*C-V(L)*S+U(L+1)
920 V(J)=U(L)*S+V(L)*C+V(L+1)
930 U(J)=U
940 NEXTJ
950 S=2*S*C:C=2*C*C-1
960 NEXTP
970 R(K)=(X(0)+(U(1)*C+V(1)*S))/N
980 NEXTK
990 REM
1000 FORK=QTOES-1
1010 J1=J2
1020 S=2*PI*K/N
1030 C=COS(S):S=SIN(S)
1040 FORJ=1TO J1
1050 L=2*J-1
1060 U(J)=-(X(L)*S)
1070 V(J)=X(L)*C+X(L+1)
1080 NEXTJ
1090 S=2*S*C:C=2*C*C-1
1100 FORP=1TOP1
1110 U(J1+1)=0:V(J1+1)=0
1120 J1=INT((J1+1)/2)
1130 FORJ=1TOJ1
1140 L=2*J-1
1150 U=U(L)*C-V(L)*S+U(L+1)
1160 V(J)=U(L)*S+V(L)*C+V(L+1)
1170 U(J)=U
1180 NEXTJ
1190 S=2*S*C:C=2*C*C-1
1200 NEXTP
1210 I(K)=-((U(1)*C+V(1)*S)/N)
1220 NEXTK
1230 REM
1240 IFES>QTHENQ=ES
1250 REM
1260 REM
1270 REM
1280 REM CALCULATE X1(0)
1290 F1=0:F2=0
1300 FORK=1TOES-1
1310 T=R(K)
1320 F1=F1+T
1330 F2=F2+K*K*T
1340 NEXTK
1350 X1(0)=R(0)+2*(F1-F2*(1/ES/ES))
1360 X1(0)=X1(0)+B
1370 REM
1380 P1=INT(100*LOG(2*ES-3)/LOG(2)+.5)/100
1390 FORJ=1TON-1
1400 T2=ES*ES
1410 FORK=1TOES-1
1420 F=1-K*K/T2
1430 U(K)=R(K)*F:V(K)=-(I(K)*F)
1440 NEXTK
1450 K1=ES-1
1460 S=2*PI*J/N
1470 C=COS(S):S=SIN(S)
1480 FORP=1TOP1
1490 U(K1+1)=0:V(K1+1)=0
1500 K1=INT((K1+1)/2)
1510 FORK=1TOK1
1520 L=2*K-1
1530 U=U(L)*C-V(L)*S+U(L+1)
1540 V(K)=U(L)*S+V(L)*C+V(L+1)
1550 U(K)=U
1560 NEXTK
1570 S=2*S*C:C=2*C*C-1
1580 NEXTP
1590 X1(J)=R(0)+2*(U(1)*C+V(1)*S)
1600 X1(J)=X1(J)+M*J+B
1610 NEXTJ
1620 RETURN